This online ressource is there to provide additional interactive visualization tool for the Varroa mites sampling and sequencing data obtained from the associated paper “XXXX”. A Snakemake pipeline created and used for genomics mapping, variant calling and analysis is provided in this host GitHub repository. Most genomics and demographic analysis using fastsimcoal2 were computed on the Okinawa Institute of Science and Technology OIST cluster, Sango.
Click on “code” button on the right to make appear the libraries used and R code for each section.
setwd("~/Dropbox (OIST)/varroa-host-switch-demo/R_Markdown")
library("tidyverse")
library("knitr") # for the markdown
library("kableExtra") # for creating a scrolling table
library("ggplot2") # for plotting
library("maps")
library("leaflet")
library("htmltools")
library("rgdal")
library("gridExtra")
library("gghighlight")
library("adegenet")
library("vcfR")
library("ggmap")
library("poppr")
library("RColorBrewer")
library("gridExtra")
library("igraph")
library("StAMPP")
library("lattice")
library("reshape2")
library("ggrepel") # for text labels
library("pcadapt")
library("data.table")
library("ggthemes")
library("ggpubr")
library("boot")
library("plotly")
library("forcats")
library("plotly")
Here, we have 43 samples of Varroa mites for which we have sequenced the whole genome succesfully. We use the leaflet package to create an interactive map showing the location and from which host, each Varroa mite was collected.
ID = Name code of the samples obtained from the CSIRO varroa collection
Species = Identified species according to mtDNA of the mtDNA COX1 partial gene
Host collected = Varroa mites were collected from within honey bee colonies identified by the collector as either Western honey bee Apis mellifera or Eastern honey bee Apis cerana.
Location/Country/Date = Exact informations on sampling
Site = Code name of each sampling site corresponding on the map Approximate.X/Y = GPS coordinates from samples obtained before 2008 are infered from approximate location given. All samples geo-coordinates after 2008 were obtained directly from field collection.
COI_haplogroup = mtDNA haplogroup identified from mtDNA COX1 partial region (reconstruct from HiSeq4000 reads and confirmed by Sanger sequencing).
Haplotype = mtDNA haplotype name by concatenating four genes (COX1, COX3, ATP6 and ND2)
Host_read = Host DNA identified by number of reads mapped against either A. cerana or A. mellifera
reads_mellifera/cerana = Number of reads mapped against reference genome of A. cerana or A. mellifera (Q > 20)
# Import the metadata for the Varroa samples
metadata <- read.csv("R_data/metadata_varroa_asia_15052020.csv", header = TRUE)
#head(metadata)
## Create a scrolling table
kable(cbind(metadata)) %>%
kable_styling(bootstrap_options = "striped", font_size = 9) %>%
scroll_box(width = "100%", height = "400px") #%>%
| ID | SEQ_ID | Species | Host_collected | Location | Country | Date | Site | Approximate.X | Approximate.Y | COI_haplogroup | Haplotype | Host_read | reads_cerana | reads_mellifera | Mean_Read_Depth | Total_reads | Mapped | Mapping_rate | Mapped_in_pair | Platform |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| V212 | VD212 | V. destructor | A. cerana | Gwangju, Jeonlabug Province | South Korea | 01/1996 | KOR02 | 35.159540 | 126.85260 | VD Korea K1 | VD Korea K1-1/K1-2 | cerana | 56 | 0 | 16.6 | 66272622 | 41762260 | 63.0 | 38388417 | HiSeq 4000 |
| V642_1 | VD642_1 | V. destructor | A. cerana | Zhuhai county, Guangdong Province | China | 11/2001 | CHN08 | 22.614010 | 112.63183 | VD China C1 | VD China C1-3 | cerana | 2857 | 7 | 47.1 | 119107269 | 116682672 | 98.0 | 112230978 | HiSeq 4000 |
| V657_1 | VD657_1 | V. destructor | A. cerana | Zhongshan, Guangdong Province | China | 04/2002 | CHN10 | 23.946090 | 114.69726 | VD China C1 | VD China C1-3 | cerana | 46 | 38 | 21.5 | 54511467 | 53363635 | 97.9 | 49849866 | HiSeq 4000 |
| V651_1 | VD651_1 | V. destructor | A. cerana | Dayao county, Yunnan Province | China | 04/2002 | CHN05 | 25.729510 | 101.33661 | VD China C3 | VD China C3-2 | cerana | 167 | 6 | 28.2 | 71261265 | 69804599 | 98.0 | 66308649 | HiSeq 4000 |
| V577_1 | VD577_1 | V. destructor | A. cerana | Yunnan Agricultural University, Kunming, Yunnan Province | China | 05/2002 | CHN06 | 25.125650 | 102.74318 | VD China C4 | VD China C4-2 | cerana | 125 | 4 | 35.0 | 91721740 | 89577965 | 97.7 | 85757545 | HiSeq 4000 |
| V654_1 | VD654_1 | V. destructor | A. cerana | Xishuangbanna, Yunnan Province | China | 04/2002 | CHN07 | 22.008810 | 100.79714 | VD Viet Nam V1 | VD Viet Nam V1-5 | cerana | 626 | 0 | 25.1 | 63453379 | 62161106 | 98.0 | 58776768 | HiSeq 4000 |
| V658_2 | VD658_2 | V. destructor | A. cerana | Xishuangbanna, Yunnan Province | China | 04/2002 | CHN07 | 22.008810 | 100.79714 | VD Viet Nam V1 | VD Viet Nam V1-6 | cerana | 2693 | 8 | 15.0 | 44339428 | 38654491 | 87.2 | 35187264 | HiSeq 4000 |
| V676_2 | VD676_2 | V. destructor | A. cerana | Chiang Mai University | Thailand | 01/2003 | THA01 | 18.805960 | 98.95336 | VD Viet Nam V1 | VD Viet Nam V1-7 | cerana | 366 | 10 | 5.9 | 29358273 | 19254785 | 65.6 | 17719549 | HiSeq 4000 |
| V149 | VD149 | V. destructor | A. cerana | Phjong, Ninh Binh Province | Viet Nam | 10/1996 | VNM01 | 21.181970 | 106.00161 | VD Viet Nam V1 | VD Viet Nam V1-8 | cerana | 1773 | 76 | 10.7 | 46623619 | 29994005 | 64.3 | 26399304 | HiSeq 4000 |
| V475_1 | VD475_1 | V. destructor | A. cerana | Hoc Chau | Viet Nam | 11/1998 | VNM02 | 20.922080 | 104.75209 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 6 | 555 | 19.3 | 50827244 | 48381793 | 95.2 | 44842843 | HiSeq 4000 |
| V153_2 | VD153_2 | V. destructor | A. mellifera | Jungup | South Korea | 10/1996 | KOR01 | 35.569880 | 126.85589 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 0 | 5824 | 22.2 | 56113114 | 55021332 | 98.1 | 51939126 | HiSeq 4000 |
| V159_1 | VD159_1 | V. destructor | A. mellifera | Jungup | South Korea | 10/1996 | KOR01 | 35.569880 | 126.85589 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 14 | 5707 | 24.2 | 64315412 | 60528716 | 94.1 | 55776173 | HiSeq 4000 |
| V639_1 | VD639_1 | V. destructor | A. mellifera | Heilongjiang Province | China | 09/2001 | CHN01 | 46.617160 | 131.15727 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 20 | 1096 | 17.2 | 45259621 | 43071591 | 95.2 | 39536470 | HiSeq 4000 |
| V625_2 | VD625_2 | V. destructor | A. mellifera | Fengman, Lishu Co., Changchun, Jilin Province | China | 05/2001 | CHN02 | 43.821600 | 126.56227 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 2 | 217 | 16.6 | 46030858 | 41845105 | 90.9 | 37822659 | HiSeq 4000 |
| V622_1 | VD622_1 | V. destructor | A. mellifera | Apiary of Jinzhong, Yuci, Shanxi Province | China | 06/2001 | CHN03 | 37.697790 | 112.70822 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 0 | 1243 | 13.8 | 39024603 | 35070629 | 89.9 | 31675460 | HiSeq 4000 |
| V641_1 | VD641_1 | V. destructor | A. mellifera | Mianzhu County, Sichuan Province | China | 09/2001 | CHN04 | 31.338070 | 104.22074 | VD Korea K1 | VD Korea K1-7 | mellifera | 6 | 871 | 21.0 | 54493356 | 52402973 | 96.2 | 48743391 | HiSeq 4000 |
| V646_1 | VD646_1 | V. destructor | A. mellifera | Zhongshan, Guangdong Province | China | 11/2001 | CHN09 | 23.927990 | 113.15883 | VD Korea K1 | VD Korea K1-7 | mellifera | 4 | 322 | 25.8 | 66768826 | 64203224 | 96.2 | 59914828 | HiSeq 4000 |
| V150_2 | VD150_2 | V. destructor | A. mellifera | Nho Quan | Viet Nam | 10/1996 | VNM03 | 20.322630 | 105.74960 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 26 | 218 | 33.2 | 86978010 | 84333771 | 97.0 | 80000889 | HiSeq 4000 |
| V474_1 | VD474_1 | V. destructor | A. mellifera | Hoc Chau | Viet Nam | 01/1998 | VNM02 | 20.922080 | 104.75209 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 52 | 648 | 24.8 | 63242091 | 61585346 | 97.4 | 58056817 | HiSeq 4000 |
| V028 | VJ028 | V. jacobsoni | A. cerana | UPM, Serdang, Kuala Lumpur | Malaysia | 03/1995 | MYS01 | 2.987890 | 101.73517 | VJ Malaysia 2 | VJ Malaysia 2-1 | cerana | 168 | 66 | 50.4 | 131769709 | 125151456 | 95.0 | 120488777 | HiSeq 4000 |
| V780 | VJ780 | V. jacobsoni | A. cerana | Chauthanh Village, Bentra District | Viet Nam | 03/2003 | VNM06 | 10.304830 | 106.35520 | VJ Laos L1 | VJ Laos L1-3 | cerana | 40372 | 4 | 15.6 | 47957642 | 38953995 | 81.2 | 36880992 | HiSeq 4000 |
| V787_2 | VJ787_2 | V. jacobsoni | A. cerana | Rach Gia Town, Kien Giang | Viet Nam | 08/2004 | VNM04 | 10.021500 | 105.09109 | VJ Laos L1 | VJ Laos L1-4 | cerana | 2325 | 28 | 29.6 | 78330554 | 75988942 | 97.0 | 72622136 | HiSeq 4000 |
| V788_1 (replicate discarded) | VJ788_1 | V. jacobsoni | A. cerana | My tho City, Kien Giang | Viet Nam | 08/2004 | VNM05 | 10.376520 | 106.34388 | VJ Laos L1 | VJ Laos L1-4 | cerana | 1291 | 2 | 7.0 | 43611088 | 19916649 | 45.7 | 18505755 | HiSeq 4000 |
| V788_2 | VJ788_2 | V. jacobsoni | A. cerana | My tho City, Kien Giang | Viet Nam | 08/2004 | VNM05 | 10.376520 | 106.34388 | VJ Laos L1 | VJ Laos L1-4 | cerana | 72 | 2 | 10.8 | 30299193 | 27206968 | 89.8 | 25941038 | HiSeq 4000 |
| V789_1 | VJ789_1 | V. jacobsoni | A. cerana | Giong Trom Town, Ben Tre | Viet Nam | 08/2004 | VNM07 | 10.157180 | 106.50445 | VJ Laos L1 | VJ Laos L1-4 | cerana | 859 | 4 | 23.6 | 71558906 | 58473472 | 81.7 | 56649683 | HiSeq 4000 |
| V333_1 | VJ333_1 | V. jacobsoni | A. cerana | Parung Panjang, Java | Indonesia | 06/1998 | IDN01 | -6.348140 | 106.56935 | VJ Java 1 | VJ Java 1-1 | cerana | 21 | 10 | 17.9 | 46019459 | 44634047 | 97.0 | 42592178 | HiSeq 4000 |
| V341_1 | VJ341_1 | V. jacobsoni | A. cerana | Sukabumi, West Java | Indonesia | 06/1998 | IDN02 | -6.923870 | 106.93156 | VJ Java 1 | VJ Java 1-2 | cerana | 119 | 4 | 24.8 | 63057170 | 61695629 | 97.8 | 59251392 | HiSeq 4000 |
| V347_1 | VJ347_1 | V. jacobsoni | A. cerana | Cililin | Indonesia | 06/1998 | IDN04 | -6.983450 | 107.43608 | VJ Java 1 | VJ Java 1-3 | cerana | 61 | 4 | 18.2 | 48579482 | 45809925 | 94.3 | 42898293 | HiSeq 4000 |
| V363_1 | VJ363_1 | V. jacobsoni | A. cerana | Johor, East Java | Indonesia | 06/1998 | IDN05 | -7.223850 | 112.73320 | VJ Java 1 | VJ Java 1-4 | cerana | 5115 | 64 | 21.3 | 54026529 | 52825280 | 97.8 | 50723792 | HiSeq 4000 |
| V565_2 | VJ565_2 | V. jacobsoni | A. cerana | Gunung Arca, Sukabumi | Indonesia | 04/2002 | IDN03 | -6.926560 | 106.95533 | VJ Java 1 | VJ Java 1-5 | cerana | 714 | 28 | 42.2 | 107959838 | 106090376 | 98.3 | 102687165 | HiSeq 4000 |
| V325_1 | VJ325_1 | V. jacobsoni | A. cerana | Kuta, Lombok | Indonesia | 04/1998 | IDN06 | -8.880270 | 116.28361 | VJ Lombok 2 | VJ Lombok 2-1 | cerana | 152 | 0 | 24.7 | 62599084 | 61346630 | 98.0 | 58922351 | HiSeq 4000 |
| V950_1 | VJ950_1 | V. jacobsoni | A. cerana | Aiyura | Papua New Guinea | 05/1997 | PNG01 | -6.343340 | 145.90809 | VJ PNG 1 | VJ PNG 1-1 | cerana | 475 | 68 | 25.8 | 66335209 | 64182154 | 96.8 | 61215275 | HiSeq 4000 |
| V854_1 | VJ854_1 | V. jacobsoni | A. cerana | Goroka | Papua New Guinea | 05/2008 | PNG04 | -6.083470 | 145.38626 | VJ PNG 1 | VJ PNG 1-1 | cerana | 5230 | 12 | 16.7 | 44613684 | 42039605 | 94.2 | 39498649 | HiSeq 4000 |
| V853_4 | VJ853_4 | V. jacobsoni | A. cerana | Daru | Papua New Guinea | 06/2008 | PNG11 | -9.078190 | 143.20997 | VJ PNG 1 | VJ PNG 1-2 | cerana | 200 | 60 | 10.4 | 35177854 | 32209884 | 91.6 | 29491225 | HiSeq 4000 |
| V983_1 | VJ983_1 | V. jacobsoni | A. cerana | Asaro | Papua New Guinea | 11/2015 | PNG05 | -6.012000 | 145.32500 | VJ PNG 1 | VJ PNG 1-1 | cerana | 284 | 102 | 19.1 | 50283170 | 47881277 | 95.2 | 44989750 | HiSeq 4000 |
| V847 | VJ847 | V. jacobsoni | A. mellifera | Wabag, Enga Province | Papua New Guinea | 06/2008 | PNG10 | -5.482670 | 143.66373 | VJ PNG 1 | VJ PNG 1-1 | mellifera | 23 | 115346 | 22.5 | 64504980 | 56218551 | 87.2 | 52755225 | HiSeq 4000 |
| V852_3 | VJ852_3 | V. jacobsoni | A. mellifera | Goroka, Eastern Highlands Province | Papua New Guinea | 05/2008 | PNG04 | -6.083470 | 145.38626 | VJ PNG 1 | VJ PNG 1-1 | mellifera | 4 | 2288 | 20.9 | 55816067 | 54523224 | 97.7 | 51390551 | HiSeq 4000 |
| V857_2 | VJ857_2 | V. jacobsoni | A. mellifera | Benabena, Eastern Highlands Province | Papua New Guinea | 05/2008 | PNG03 | -6.128125 | 145.42904 | VJ PNG 1 | VJ PNG 1-1 | mellifera | 20 | 429 | 18.8 | 57477057 | 46759609 | 81.4 | 44551404 | HiSeq 4000 |
| V856_1 | VJ856_1 | V. jacobsoni | A. mellifera | Henganofi Border, Eastern Highlands Province | Papua New Guinea | 05/2008 | PNG02 | -6.260350 | 145.59681 | VJ PNG 1 | VJ PNG 1-1 | mellifera | 15 | 9803 | 13.6 | 58525652 | 37045620 | 63.3 | 31803798 | HiSeq 4000 |
| V952_1 | VJ952_1 | V. jacobsoni | A. mellifera | Amayufa, Eastern Highlands Province | Papua New Guinea | 05/2014 | PNG06 | -5.975581 | 145.27676 | VJ PNG 1 | VJ PNG 1-1 | mellifera | 18 | 36 | 22.3 | 61969612 | 55519034 | 89.6 | 52870625 | HiSeq 4000 |
| V953_2 | VJ953_2 | V. jacobsoni | A. mellifera | Benabena, Eastern Highlands Province | Papua New Guinea | 05/2014 | PNG03 | -6.128125 | 145.42904 | VJ PNG 1 | VJ PNG 1-1 | mellifera | 22 | 26196 | 22.3 | 62054570 | 56103596 | 90.4 | 53032044 | HiSeq 4000 |
| V954_2 | VJ954_2 | V. jacobsoni | A. mellifera | Windy Lodge | Papua New Guinea | 06/2014 | PNG07 | -6.035146 | 144.96032 | VJ PNG 1 | VJ PNG 1-1 | mellifera | 14 | 159 | 22.5 | 58553267 | 56180797 | 96.0 | 52846232 | HiSeq 4000 |
| V955_1 | VJ955_1 | V. jacobsoni | A. mellifera | Donakanage | Papua New Guinea | 06/2014 | PNG08 | -5.882240 | 144.82077 | VJ PNG 1 | VJ PNG 1-1 | mellifera | 15 | 460 | 24.7 | 64865555 | 61577522 | 94.9 | 58553870 | HiSeq 4000 |
| V956_4 | VJ956_4 | V. jacobsoni | A. mellifera | Minj | Papua New Guinea | 06/2014 | PNG09 | -5.826640 | 144.40741 | VJ PNG 1 | VJ PNG 1-2 | cerana | 1444 | 69 | 17.1 | 47824500 | 45966716 | 96.1 | 42702569 | HiSeq 4000 |
| Mom #1 | Mom #1 | V. destructor | A. mellifera | OIST Ecology and Evolution experimental apiary, Okinawa | Japan | 02/2018 | JPN01 | 26.456448 | 127.82386 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 0 | 7448 | 20.2 | 52944228 | 50821794 | 96.0 | 52944228 | HiSeq 4000 |
| Son #1 | Son #1 | V. destructor | A. mellifera | OIST Ecology and Evolution experimental apiary, Okinawa | Japan | 02/2018 | JPN01 | 26.456448 | 127.82386 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 0 | 341 | 19.8 | 51620028 | 49266370 | 95.4 | 51620028 | HiSeq 4000 |
| Mom #2 | Mom #2 | V. destructor | A. mellifera | OIST Ecology and Evolution experimental apiary, Okinawa | Japan | 02/2018 | JPN01 | 26.456448 | 127.82386 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 0 | 50765 | 24.8 | 65220912 | 62435130 | 95.7 | 65220912 | HiSeq 4000 |
| Son #2 | Son #2 | V. destructor | A. mellifera | OIST Ecology and Evolution experimental apiary, Okinawa | Japan | 02/2018 | JPN01 | 26.456448 | 127.82386 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 0 | 422 | 18.7 | 48499085 | 46511492 | 95.9 | 48499085 | HiSeq 4000 |
| Mom #7 | Mom #7 | V. destructor | A. mellifera | OIST Ecology and Evolution experimental apiary, Okinawa | Japan | 02/2018 | JPN01 | 26.456448 | 127.82386 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 0 | 563 | 16.3 | 43293860 | 41401282 | 95.6 | 43293860 | HiSeq 4000 |
| Son #7 | Son #7 | V. destructor | A. mellifera | OIST Ecology and Evolution experimental apiary, Okinawa | Japan | 02/2018 | JPN01 | 26.456448 | 127.82386 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 0 | 660 | 18.8 | 49134877 | 46899128 | 95.5 | 49134877 | HiSeq 4000 |
| Mom #15 | Mom #15 | V. destructor | A. mellifera | OIST Ecology and Evolution experimental apiary, Okinawa | Japan | 05/2018 | JPN01 | 26.456448 | 127.82386 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 0 | 100 | 22.5 | 59835146 | 56138955 | 93.8 | 59835146 | HiSeq 4000 |
| Son #15 | Son #15 | V. destructor | A. mellifera | OIST Ecology and Evolution experimental apiary, Okinawa | Japan | 05/2018 | JPN01 | 26.456448 | 127.82386 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 1 | 193 | 21.4 | 57025270 | 53922121 | 94.6 | 57025270 | HiSeq 4000 |
| Mom #17 | Mom #17 | V. destructor | A. mellifera | OIST Ecology and Evolution experimental apiary, Okinawa | Japan | 05/2018 | JPN01 | 26.456448 | 127.82386 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 2 | 40283 | 22.1 | 58014781 | 55352009 | 95.4 | 58014781 | HiSeq 4000 |
| Son #17 | Son #17 | V. destructor | A. mellifera | OIST Ecology and Evolution experimental apiary, Okinawa | Japan | 05/2018 | JPN01 | 26.456448 | 127.82386 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 2 | 949 | 33.0 | 87139528 | 82953771 | 95.2 | 87139528 | HiSeq 4000 |
| Mom #19 | Mom #19 | V. destructor | A. mellifera | OIST Ecology and Evolution experimental apiary, Okinawa | Japan | 05/2018 | JPN01 | 26.456448 | 127.82386 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 0 | 24181 | 18.1 | 47748745 | 45343415 | 95.0 | 47748745 | HiSeq 4000 |
| Son #19 | Son #19 | V. destructor | A. mellifera | OIST Ecology and Evolution experimental apiary, Okinawa | Japan | 05/2018 | JPN01 | 26.456448 | 127.82386 | VD Korea K1 | VD Korea K1-1/K1-2 | mellifera | 5 | 498 | 24.5 | 65304925 | 61543197 | 94.2 | 65304925 | HiSeq 4000 |
#row_spec(1:23, bold = T, color = "white", background = "#E60000") %>%
#row_spec(24:39, bold = T, color = "white", background = "#FF9999") %>%
#row_spec(40:74, bold = T, color = "white", background = "#0000FF") %>%
#row_spec(75:96, bold = T, color = "white", background = "#3399FF")
# Download the country borders layer
#download.file("http://thematicmapping.org/downloads/TM_WORLD_BORDERS_SIMPL-0.3.zip" , destfile="world_shape_file.zip")
#system("unzip world_shape_file.zip")
world_spdf=readOGR(dsn= getwd() , layer="TM_WORLD_BORDERS_SIMPL-0.3")
## OGR data source with driver: ESRI Shapefile
## Source: "/Users/alphamanae/Dropbox (OIST)/varroa-host-switch-demo/R_Markdown", layer: "TM_WORLD_BORDERS_SIMPL-0.3"
## with 246 features
## It has 11 fields
## Integer64 fields read as strings: POP2005
data.vdac <- metadata %>% filter(Species == "V. destructor") %>% filter(Host_collected == "A. cerana")
data.vdam <- metadata %>% filter(Species == "V. destructor") %>% filter(Host_collected == "A. mellifera")
data.vjac <- metadata %>% filter(Species == "V. jacobsoni") %>% filter(Host_collected == "A. cerana")
data.vjam <- metadata %>% filter(Species == "V. jacobsoni") %>% filter(Host_collected == "A. mellifera")
The following interactive map showed the distribution of both V. destructor and V. jacobsoni samples on either their original host the eastern honey bee A. cerana or on the new host the western honey bee A. mellifera.
If you pass your mouse over a particular point, a pop-up will appear presenting essential information about the particular Varroa mite sample. On the right top corner, a layer control button allow to subselect species/host couple.
V. destructor on A. cerana
V. destructor on A. mellifera
V. jacobsoni on A. cerana
V. jacobsoni on A. mellifera
## Create the interactive map with leaflet
leaflet(metadata) %>%
addTiles(group = "OSM (default)") %>%
## add the layer other than default we would like to use for background
addProviderTiles(providers$CartoDB.PositronNoLabels, group = "Positron NoLabels") %>%
## adding the layers with V. destructor
addCircleMarkers(data = data.vdac, data.vdac$Approximate.Y, data.vdac$Approximate.X,
weight = 0.5,
col = "#FB0000",
radius = 4,
fillOpacity = 0.9,
stroke = T,
popup = ~paste(sep = "<br/>",
paste ("Name: ", as.character(ID)),
paste ("Host: ", as.character(Host_collected)),
paste ("Date: ", as.character(Date)),
paste ("mtDNA: ", as.character(Haplotype))),
group = 'V. destructor on A. cerana') %>%
addCircleMarkers(data = data.vdam, data.vdam$Approximate.Y, data.vdam$Approximate.X,
weight = 0.5,
col = "#FF9999",
radius = 4,
fillOpacity = 0.9,
stroke = T,
popup = ~paste(sep = "<br/>",
paste ("Name: ", as.character(ID)),
paste ("Host: ", as.character(Host_collected)),
paste ("Date: ", as.character(Date)),
paste ("mtDNA: ", as.character(Haplotype))),
group = 'V. destructor on A. mellifera') %>%
## adding the layers with V. jacobsoni
addCircleMarkers(data = data.vjac, data.vjac$Approximate.Y, data.vjac$Approximate.X,
weight = 0.5,
col = "#0000FF",
radius = 4,
fillOpacity = 0.9,
stroke = T,
popup = ~paste(sep = "<br/>",
paste ("Name: ", as.character(ID)),
paste ("Host: ", as.character(Host_collected)),
paste ("Date: ", as.character(Date)),
paste ("mtDNA: ", as.character(Haplotype))),
group = 'V. jacobsoni on A. cerana') %>%
addCircleMarkers(data = data.vjam, data.vjam$Approximate.Y, data.vjam$Approximate.X,
weight = 0.5,
col = "#00CCFF",
radius = 4,
fillOpacity = 0.9,
stroke = T,
popup = ~paste(sep = "<br/>",
paste ("Name: ", as.character(ID)),
paste ("Host: ", as.character(Host_collected)),
paste ("Date: ", as.character(Date)),
paste ("mtDNA: ", as.character(Haplotype))),
group = 'V. jacobsoni on A. mellifera') %>%
## adding the control button to remove or add layers of points
addLayersControl(position = "topright",
baseGroups = c("OSM (default)", "Positron NoLabels"),
overlayGroups = c("V. destructor on A. cerana",
"V. destructor on A. mellifera",
"V. jacobsoni on A. cerana",
"V. jacobsoni on A. mellifera"),
options = layersControlOptions(collapsed = TRUE)) %>%
## show the positron background prerably to the OSM layer
showGroup("Positron NoLabels")
Here we plot the mean read depth computed from each .bam file mapped to the reference vdes3.0.fasta using samtools depth -a FILE.bam | awk '{c++;s+=$3}END{print s/c}' in relation to the mapping percentage.
As you can see by moving your pointer on the graph, the individual V788_1 has been removed from the analysis and replaced by the replicate V788_2 with much better reads statistics.
raichu <- ggplot(metadata, aes(x=Mapping_rate, y = Mean_Read_Depth, col=Species, shape=Host_collected, label = ID))
mycol = c("red2", "blue2")
raichu <- raichu + geom_point(size = 3)
raichu <- raichu + scale_color_manual(values = mycol)
raichu <- raichu + scale_size_manual(values = mycol)
raichu <- raichu + theme_calc()
raichu <- raichu + xlab("Mapping rate to Vdes_3.0 reference genome(%)")
raichu <- raichu + ylab("Mean average read depth")
raichu <- raichu + ylim(0,20)
#raichu
## make an interactive version of the scatter plot
intraichu <- ggplotly(raichu)
intraichu
We highlight the mean read depth regarding the general distribution, where is each group of individuals “species per host” located in it.
evoli <- ggplot(metadata, aes(Mean_Read_Depth, fill = Species))
evoli <- evoli + geom_histogram(bins = 44)
evoli <- evoli + gghighlight()
evoli <- evoli + facet_wrap(~Species~Host_collected)
evoli <- evoli + scale_fill_manual(values = c("red", "blue"))
evoli <- evoli + theme_classic()
evoli <- evoli + labs(x = "Mean Read Depth using NextGenMap", y = "Number of samples")
evoli <- evoli + ggtitle("Read depth does not appear as unbalanced accross species and host")
evoli
Following is a plot of the total number of reads generated per individual and the proportion of reads mapped onto the reference.
togepi <- ggplot(metadata, aes(x=ID, label = Mapping_rate))
togepi <- togepi + geom_bar(aes(weight=Total_reads), fill="black", position="dodge")
togepi <- togepi + geom_bar(aes(weight=Mapped), fill="#1d96f2", position="dodge")
togepi <- togepi + theme(axis.text.x = element_text(angle = 90, hjust = 1))
togepi <- togepi + labs(x = "Individuals", y = "Number of reads")
togepi <- togepi + ggtitle("Proportion of reads mapped against reference using NextGenMap")
togepi
Before proceeding to DNA extraction with a destructive method, we measured the body size of each Varroa mite whenever possible on both dorsal and ventral view. We measured seven morphological characters but retained only the body width and length for comparison here between Varroa species and their associated honey bee hosts.
The scatterplot also shows the average values obtained by Anderson and Trueman (2000) Exp. App. Acar. for V. destructor, V. jacobsoni and undetermined mites from Luzon and Mindanao.
# Table with 512 morpho data for the moment
morpho <- read.csv("R_data/Varroa-morpho-subset.csv", header = TRUE)
#kable(cbind(morpho)) %>%
# kable_styling(bootstrap_options = "striped", font_size = 10) %>%
# scroll_box(width = "100%", height = "400px")
morpho %>%
plot_ly(x = ~BW_V, y = ~BL_V,
color = ~Species,
symbol = ~Host,
hoverinfo = "text",
text = ~paste("Name: ", Sample.ID, "<br>",
"Species: ", Species, "<br>",
"Host: ", Host, "<br>",
"Country: ", Country.Island)) %>%
add_markers(colors = c("indianred2", "red3", "skyblue2", "blue2", "black"),
symbols = c( "triangle-up", "circle", "x"),
size = 3) %>%
layout(xaxis = list(title = "Body width (mm)"),
yaxis = list(title = "Body length (mm)"),
title = "Variation in body size measured from ventral view")
morpho %>%
plot_ly(x = ~BW_D, y = ~BL_D,
color = ~Species,
symbol = ~Host,
hoverinfo = "text",
text = ~paste("Name: ", Sample.ID, "<br>",
"Species: ", Species, "<br>",
"Host: ", Host, "<br>",
"Country: ", Country.Island)) %>%
add_markers(colors = c("indianred2", "red3", "skyblue2", "blue2", "black"),
symbols = c( "triangle-up", "circle", "x"),
size = 3) %>%
layout(xaxis = list(title = "Body width (mm)"),
yaxis = list(title = "Body length (mm)"),
title = "Variation in body size measured from dorsal view")
STATS TEST?
After formatting the variant call file .vcf into a plink format by renaming each chromosome (e.g., NW_019211454.1 by 1, NW_019211455.1 by 2, …), and applying the option pca, we obtained the two files .eigenval and .eigenvec.
The following 2D and 3D PCAs are based on the ~1.4 million biallelic SNPs with a minor allelic frequency higher than 5%.
pca.all<- read_table2("R_data/43indplink.eigenvec", col_names = FALSE)
eigenval.all <- scan("R_data/43indplink.eigenval")
## Here the first two columns contains the individuals names as we did not specify it in plink earlier
# We remove the first column and assign the new name of column 1
pca.all <- pca.all[,-1]
names(pca.all)[1] <- "ID_NAMES"
# Rename the columns for PCA axis
names(pca.all)[2:ncol(pca.all)] <- paste0("PC", 1:(ncol(pca.all)-1))
## Joint both table by using the FILE_ID as common parameters
metapca.all <- left_join(pca.all, metadata, by = c("ID_NAMES" = "SEQ_ID"))
### PLOTS
# first convert to percentage variance explained
pve <- data.frame(PC = 1:20, pve = eigenval.all/sum(eigenval.all)*100)
# make Eigenplot
a <- ggplot(pve, aes(PC, pve)) + geom_bar(stat = "identity")
a + ylab("Percentage variance explained") + theme_light()
metapca.all %>%
plot_ly(x = ~PC1, y = ~PC2,
color = ~Species,
symbol = ~Host_collected,
hoverinfo = "text",
text = ~paste("Name: ", ID, "<br>",
"Species: ", Species, "<br>",
"Host: ", Host_collected, "<br>",
"Country: ", Country)) %>%
add_markers(colors = c("red3", "blue2"),
symbols = c( "triangle-up", "circle"),
size = 4)%>%
layout(title = "Varroa mites genetically differentiate between hosts in both species")
metapca.all %>%
plot_ly(x = ~PC1, y = ~PC2, z = ~PC3,
color = ~Country,
hoverinfo = "text",
text = ~paste("Name: ", ID, "<br>",
"Species: ", Species, "<br>",
"Host: ", Host_collected, "<br>",
"Country: ", Country,
"Haplogroup: ", COI_haplogroup)) %>%
add_markers(colors = c("red3", "purple", "blue2", "green", "orange", "brown", "pink"),
size = 10) %>%
layout(title = "Mites on A. cerana genetically differs within the native range")
metapca.all %>%
plot_ly(x = ~PC1, y = ~PC2, z = ~PC3,
color = ~COI_haplogroup,
hoverinfo = "text",
text = ~paste("Name: ", ID, "<br>",
"Species: ", Species, "<br>",
"Host: ", Host_collected, "<br>",
"Country: ", Country,
"Haplogroup: ", COI_haplogroup)) %>%
add_markers(colors = c("red3", "purple", "blue2", "green", "orange", "brown", "pink"),
size = 10) %>%
layout(title = "Mitochondrial lineages can be a pre-indicator of whole genome differentiation")
Same idea as the previous PCA, except that only 236,713 SNPs were kept to reduce the effect of linkage disequilibrium.
We can see that even with a smaller reduce number of SNPs, Varroa mites species are genetically distinct and
pca.allLD<- read_table2("R_data/43ind-ldpruned.eigenvec", col_names = FALSE)
eigenval.allLD <- scan("R_data/43ind-ldpruned.eigenval")
## Here the first two columns contains the individuals names as we did not specify it in plink earlier
# We remove the first column and assign the new name of column 1
pca.allLD <- pca.allLD[,-1]
names(pca.allLD)[1] <- "ID_NAMES"
# Rename the columns for PCA axis
names(pca.allLD)[2:ncol(pca.allLD)] <- paste0("PC", 1:(ncol(pca.allLD)-1))
## Joint both table by using the FILE_ID as common parameters
metapca.allLD <- left_join(pca.allLD, metadata, by = c("ID_NAMES" = "SEQ_ID"))
### PLOTS
# first convert to percentage variance explained
pve <- data.frame(PC = 1:20, pve = eigenval.allLD/sum(eigenval.allLD)*100)
# make Eigenplot
a <- ggplot(pve, aes(PC, pve)) + geom_bar(stat = "identity")
a + ylab("Percentage variance explained") + theme_light()
metapca.allLD %>%
plot_ly(x = ~PC1, y = ~PC2,
color = ~Species,
symbol = ~Host_collected,
hoverinfo = "text",
text = ~paste("Name: ", ID, "<br>",
"Species: ", Species, "<br>",
"Host: ", Host_collected, "<br>",
"Country: ", Country)) %>%
add_markers(colors = c("red3", "blue2"),
symbols = c( "triangle-up", "circle"),
size = 4) %>%
layout(title = "~200,000 LD pruned SNPs still allowed to see host genetic divergence")
metapca.allLD %>%
plot_ly(x = ~PC1, y = ~PC2, z = ~PC3,
color = ~Country,
hoverinfo = "text",
text = ~paste("Name: ", ID, "<br>",
"Species: ", Species, "<br>",
"Host: ", Host_collected, "<br>",
"Country: ", Country,
"Haplogroup: ", COI_haplogroup)) %>%
add_markers(colors = c("red3", "purple", "blue2", "green", "orange", "brown", "pink"),
size = 10) %>%
layout(title = "3D PCA color-coded by country")
metapca.allLD %>%
plot_ly(x = ~PC1, y = ~PC2, z = ~PC3,
color = ~COI_haplogroup,
hoverinfo = "text",
text = ~paste("Name: ", ID, "<br>",
"Species: ", Species, "<br>",
"Host: ", Host_collected, "<br>",
"Country: ", Country,
"Haplogroup: ", COI_haplogroup)) %>%
add_markers(colors = c("red3", "purple", "blue2", "green", "orange", "brown", "pink"),
size = 10) %>%
layout(title = "3D PCA color-coded by mtDNA COX1 haplogroup")